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Motivated by experiments on sheared suspensions that show a transition between ordered and 
disordered phases, we here study the long-time behavior of a sheared and overdamped 2-d system 
of particles interacting by repulsive forces. As a function of interaction strength and shear rate 
we find transitions between phases with vanishing and large single-particle diffusion. In the phases 
with vanishing single-particle diffusion, the system evolves towards regular lattices, usually on very 
slow time scales. Different lattices can be approached, depending on interaction strength and forcing 
amplitude. The disordered state appears in parameter regions where the regular lattices are unstable. 
Correlation functions between the particles reveal the formation of shear bands. In contrast to single 
particle densities, the spatially resolved two-particle correlation functions vary with time and allow 
to determine the phase within a period. As in the case of the suspensions, motion in the state with 
low diffusivity is essentially reversible, whereas in the state with strong diffusion it is not. 

PACS numbers: 45.70.-n, 45.50.-j, 05.45.-a, 05.65.-l-b 


I. INTRODUCTION 

The observation of echoes, i.e. a recovery of the initial 
state, for spins processing in a magnetic field upon rever¬ 
sal of the field [1, 2] and for dye in viscous fluids [3] are 
among the best known experiments that seemingly defy 
the expected irreversibility in the evolution of many body 
systems. Several other examples have followed [4, 5], and 
in the context of chaotic dynamical systems the absence 
of echoes and the inability to return to initial conditions 
has frequently been used as a test for chaotic dynam¬ 
ics [6-9]. While in all these cases the strength of the 
echoes deteriorates gradually as parameters are varied, 
the experiments by Pine et al [10] on sheared suspen¬ 
sions added a new facet to the old problem: they find a 
phase-transition like change from an essentially reversible 
state for one set of parameters to an irreversible state for 
other parameters. Subsequent studies have demonstrated 
similar transitional behavior in other hydrodynamic ex¬ 
periments [11-13], particle systems [14, 15], and comple¬ 
mentary simulations [16] and simplified models [17-21]. 
It has also been put into the context of more general 
order-disorder transitions in solid materials [22-24] and 
granular systems [25-27]. Similar transitions were also 
reported for superconducting vortices [28-30]. 

The demonstration of the transition in several other 
systems suggests that it is a general phenomenon for 
sheared many body systems. This is also supported by 
the ability to capture much of the observed behavior and 
the transition in the discrete time model described in [18]. 
In order to explore the origin of the phenomenon further 
and to also establish a connection to reversibility studies 
in chaotic dynamical systems, we here consider a model 
for a many-body system that contains some of the fea¬ 
tures of the suspension experiment but allows for a much 
more detailed investigation of its long-time properties. 

In section II we present the details of the model and 
demonstrate that it also shows the transition observed 


in [10]. In section III we characterize the system both 
before and after the disorder transition. In section IV 
we investigate the stability properties of the limit cycles 
described in the previous sections. Finally, in section V 
we analyze the spatially resolved correlations and their 
phase dependency. We conclude with a few general re¬ 
marks and observations in section VI. 

II. THE MODEL AND ITS PARAMETERS 

For our model we do not aim to calculate the de¬ 
tailed hydrodynamic interaction or exact forces between 
charged colloids, as done in the simulations of [10]. 
Rather, we pick an interaction force that captures key 
features of the fluid-particle system. This is similar in 
spirit to the model in [18], but differs in that the model 
is a continuous one and not an event-driven mapping. 
We assume that the motion of the many body system is 
overdamped, so that there is no motion if the external 
forcing ceases. The forces on the particles are hence bal¬ 
anced by viscous friction, and the equations of motion 
become 


Here, /r is a friction coefficient, the position of parti¬ 
cle z, and Fi are the forces acting on the particle. We 
take the interaction between the particles as repulsive, 
thereby mimicking the repulsion due to liquid pressure 
when two particles come close. The particles are confined 
to a plane, and the domain is taken to be rectangular, and 
not curved as between the cylinders in the experiments of 
Pine et al. [10]. The domain is periodically continued in 
both directions. The main control parameter will be the 
amplitude of the forcing. A second parameter controls 
the strength of the interaction between the particles; it 
is also related to the density of the system, as will be 
discussed below. 
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A. Forces 


The repulsion between the particles is modeled as a po¬ 
tential with power-law decay, V(rij) oc with the 
distance between particles i and j. We here take a = 2 
as a compromise between strong repulsion for short dis¬ 
tances (i.e. larger a) and a numerically controllable range 
of the interaction (where hydrodynamics would suggest a 
decay as slow as a = 1). With this repulsive interaction, 
the system is reminiscent of Wigner crystals and clusters 
of charged particles in plasmas [31]. Particle positions 
are denoted = (xi,yi), with x along the direction of 
shear and y normal to it. The total force acting on the 
i-th particle is the sum of the mutual interactions with all 
other particles in the system and a periodic shear force in 
the x-direction whose amplitude increases linearly with 
Vi, 


F, = -V 




SyiCOs{ujt)ej;. ( 2 ) 


The parameter A determines the strength of the inter 
particle potential, Syi is the amplitude of the shear force 
and ui is its frequency. Substituted into equation (1), we 
obtain the full expression for the time evolution of the 
position of the i-th particle. 


dxj 

dt 





-k 


—j/j cos{ujt)e^ 


( 3 ) 


B. Parameters 


In order to expose the independent parameters of the 
system, we introduce a length scale A and a time scale T, 
and set x^ = Ax' and t = Tt'. The time T can naturally 
be identified with the period of the shear, T = 27r/a; . 
The length scale is not related to a quantity explicitly 
displayed in the equations, but it enters implicitly in the 
many-body system via the mean distance between parti¬ 
cles, and is thus related to the density. Therefore, vari¬ 
ations in period and density are absorbed into the two 
remaining parameters of the system, the dimensionless 
interaction strength 


Ti = 2ATly.\^ (4) 

and the dimensionless shear rate 


Ts = STjy.. 

The evolution equations then become 


dx] 

dt 


riE 



Tsy'iCos{2TTt)e^. 


( 5 ) 

( 6 ) 


We will drop the primes in the subsequent sections and 
refer to this equation as the evolution equation. 


The instantaneous strain 7(t) follows from the time 
evolution of the distance between two particles that ini¬ 
tially are displaced by Ay perpendicular to the shear. 
The separation in x-direction is then given by Ax(t) = 
'y{t)Ay with 

lit) = = rs/(27r) sin(27rt) = 70 sin(27rt); (7) 

Ay 

7o = rs/(27r) is the amplitude of the affine shear. This 
also gives a notion of the shear process: At the start of 
the period, t = 0, the strain is zero. Over the next half 
period, the system is tilted to the right, reaching its apex 
at t = r/4. Here, the flow comes to a halt and reverses 
its direction. After half a period, strain is again zero and 
the process repeats to the left. The total accumulated 
strain over one period then becomes 

7 = [ iTs cos{27rt)\ dt = 470. (8) 

C. Numerical implementation 

In order to solve equation (6), we introduce Lees- 
Edwards boundary conditions [32] to account for the 
sheared images in the y-direction. We use a modified 
minimum image convention, taking several closest images 
of each particle into account, typically one or two in each 
direction. The typical number of particles in the main 
domain is about 100, and in a few examples we went up 
to 900 particles. Time integration is performed using a 
4th order Runge-Kutta-Fehlberg integrator provided by 
the GNU Scientific Library (GSL). Computation of the 
right-hand side is parallelized using openMP. The width 
(in the x-direction) of the box is K, whereas the height 
(in the y-direction) of the box is sin(7r/3)L with AT, L € N 
and even and KL = N, as suggested by a regular triangu¬ 
lar lattice. The simulation was carried out with random 
initial conditions, and the motion was typically followed 
over several thousand periods in order to obtain clear ev¬ 
idence for the asymptotic state and to extract reliable 
statistics in the case of chaotic states. 


D. Reduced representations and stroboscopic maps 

Even for small shear, particles undergo large scale mo¬ 
tions along the shear axis, with an amplitude that in¬ 
creases with the displacement in the normal direction. 
In order to remove this affine deformation, we introduce 
a reduced representation where the translation by the 
shear is taken out, i.e. we study 

Xi{t) = Xi{t) - 7oy(to) sin(27rt). (9) 

Figure 1 shows examples of such reduced motions: for the 
chosen parameters the amplitudes of the motion around 
the reference positions are small and most particles stay 
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FIG. 1. (Color online) Particle trajectories over one period. 
Initial positions are marked by open circles, final positions by 
filled red circles. Shown are the reduced positions (9) after 
subtracting the affine transformation. Non-reversible motions 
(in particular in the lower left corner) are connected with close 
encounters and large displacements, see also supplementary 
movie online [.l.l]. The parameters are Pi = 0.1 and Fa = 16.0. 


close to their initial positions, but as the group in the 
lower left corner shows, some of them can experience 
large displacements as a consequence of close encounters. 

The long time behavior of the system shows up after 
a large number of periods only, and is then best studied 
in the form of stroboscopic maps, with positions taken 
at multiples of the period, Xi(n) = Xi{nT). In figure 2 
and the accompanying movies [33] we show two exam¬ 
ples, one for smaller shear rate and one for larger shear 
rate. For moderate shear rates, only small displacements 
are observed and the particles evolve slowly towards reg¬ 
ular lattices (section III A). For larger shear, almost all 
particles travel larger distances from their initial posi¬ 
tion. This corresponds to the disordered state of Pine 
et al. [10]. The random motion of the particles can be 
captured by diffusion coefficients, which are anisotropic 
and differ in the longitudinal and transverse directions 
(section IIIB). 

III. ORDERED AND DISORDERED STATES 
A. Ordered states 

In the absence of external forcing, the interactions be¬ 
tween the particles push a random initial condition to¬ 
wards a force equilibrium. A regular lattice, such as a 
hexagon, is an example of such an equilibrium, but it 
will rarely be reached from a random initial condition 
since particles will be trapped in a state with a few dis¬ 
locations or in a state showing several oriented patches, 
separated by grain boundaries. 

As soon as shear is added to the system, the particles 
will start to rearrange and to self-organize. After suffi¬ 
ciently long times, states like the ones shown in Fig. 3 for 


(a) 




x{nT) 


FIG. 2. (Golor online) Stroboscopic map of particle trajec¬ 
tories for low shear (Fg = 8, top) and large shear (Fg = 16, 
bottom). Open circles indicate initial positions, full red cir¬ 
cles final positions. In the top diagram for small shear one 
notes that the particles move very little over the 500 time 
steps. There is a localized region with larger motions in the 
lower right corner that is connected with a rare local rear¬ 
rangement. In the lower diagram for larger shear almost all 
particles experience large displacements over a much shorter 
interval of 25 time units. In both cases Fi = 0.1. 

Fg = 2, 6, and 8 are obtained. For very low shear rates 
the system settles into the hexagonal configuration, with 
an orientation relative to the shear direction that de¬ 
pends on the initial conditions. Since the regular hexag¬ 
onal state is a possible force-equilibrium, the external 
forcing can be seen as a minor perturbation that hardly 
affects particle interactions. With increasing shear rate, 
the lattice orientation may change, and the lattice will 
align parallel or perpendicular to the external force. At 
a shear rate around 6, the asymptotic state is given by 
a rectangular lattice configuration. At Fg « 8, we again 
observe a hexagonal lattice, though this time it is almost 
always oriented parallel to the shearing motion. Increas¬ 
ing the rate further the system will at some point fail 
to approach a regular lattice, and a chaotic, disordered 
state is attained. 

Since the asymptotic reversible states are also lattice 
configurations, the external forcing only causes affine de¬ 
formations of the system. This is in contrast to the 
observations by Keim and Arratia [15] who found clus- 
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FIG. 3. (Color online) Asymptotic states at moderate shear 
rates, Fa = 2, 6 and 8, respectively. Except for minor defects 
they are either a hexagonal or a rectangular lattice configu¬ 
ration. Fi = 0.1 in all three cases. 
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ters of particles subject to non-affine deformations in the 
reversible regime. We only observe non-affine deforma¬ 
tions in the irreversible regime and in the trajectories 
approaching the ordered states. 

The regular states observed here are a consequence of 
the homogeneity of the system, the density, and the do¬ 
main size: Imperfections will most likely trap the system 
in a glassy or disordered state, as in the experiments of 
Pine et al. [10]. Another important aspect is the dif¬ 
ferent kind of interactions considered in [10, 18]. Since 
there an interaction between particles takes place upon 
contact only, reversibility is attained as soon as each par¬ 
ticle has cleared a sufficiently large space around it. In 
our case, particles interact over longer distances so that 
further ordering is required and eventually ordered states 
are attained. Nevertheless, as Corte et al. [18] pointed 
out, very similar absorbing states exist in their systems 
as well, but are not found by the simulations. Similar 
to our case, they become unstable as soon as a critical 
strain is exceeded. 

In order to quantify the ordering we investigate the 
orientation number of an ensemble, defined by 

4'fc = (10) 

Here, k is the number of next neighbors taken into ac¬ 
count, and Qij is the angle of the vector between particles 
i and j towards the x-axis. Of particular interest are the 
values of I'I'gl, which becomes unity for a hexagonal lat¬ 
tice, and of 14141 , which becomes unity for a rectangular 
lattice. 

We show the results in Fig. 4 for a wide range in param¬ 
eter space Fi x Fg. For low shear rates and independent 
of the interaction strength, the orientation number re¬ 
flects the previously mentioned change of the asymptotic 
states. For very low shear rates up to Fg ;< 3 and for shear 
rates larger than Fg > 7, a hexagonal lattice is preferred. 
For shear rates in between we find |4'4| « 1, and thus 
a rectangular lattice is attained. With increasing shear 
rate and low enough interaction strength Fj, order is lost 
beyond Fg > 11. This corresponds to the observation of 
chaotic states such as in Fig. 2(b). For a larger inter¬ 
action strength Fj, ordered states are attained even at 


FIG. 4. (Color online) Order parameters |4'4| and jtfej for 
different values of Fi and Fg. For a low shear rate Fg we 
clearly observe a change from sixfold orientation to fourfold 
orientation (at Fg « 3) and back to sixfold orientation (at 
Fg « 7). In the light area in the upper left, no ordering of any 
kind is observable. Its lower boundary is located near Fg « 11. 
For larger Fi, the system attains lattice configurations even at 
larger shear rates. Although the system then does not reach 
perfect lattice configurations, a modulation with Fg is still 
discernible. The black solid lines indicate the loss of ordering 
where both |4'4| and j’Fgj drop below 0.4. 

large shear rates. However, the computational demand 
increases and integration times become too short, and 
hence several defects remain in the configurations and 
reduce the orientation number. However, in contrast to 
the orderless region, we find an orientation number which 
is larger than 0.4 in this region. Additionally, we observe 
a modulation of the orientation number with Fg, in which 
both lattice types alternately dominate. 

B. Diffusion constants 

Convenient scalar quantities that can be used to distin¬ 
guish the different long-time dynamics are related to the 
diffusive motions of the particles. The local and instan¬ 
taneous quantity is the mean displacement of particles 
over one period, 

1 ^ 

' fc=i 

It provides a convenient measure to quickly distinguish 
between regular and irregular states, as shown in Fig. 5. 
For large shear (Fg = 16) the displacement is large and 
varies little in time, indicating a disordered, chaotic state. 
For low shear (Fg = 8) it decreases, eventually reaching 
machine precision: the system evolves towards a stable 
state. The intermediate maxima, such as the one near 
t = 2800 in Fig. 5, indicate reorganizations to remove 
lattice defects. The event is related to the rearrangement 
in the lower right corner in Fig. 2 (top) and clearly visible 
in the movie [33]. 
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FIG. 5. (Color online) Average displacement after one pe¬ 
riod for different shear rates Fg. For low shear rates, it drops 
to numerical precision with intermediate bumps indicating 
reorganizations. For large shear rates, the average displace¬ 
ment approaches a finite value reflecting the random motion. 
Fi = 0.1. 


In the diffusive state, we can calculate diffusion con¬ 
stants. In view of the asymmetry in the system in the 
directions along the shear and perpendicular to it, we 
consider two separate diffusions constants. The diffusion 
constants can be obtained either by averaging in time 
over the displacements over one period, 

2D,r7= /^^|xfc(n)-a;fc(n-l)|A (12) 

' 1 ' ^^ritransient 

or by approximating the mean square displacement of the 
particles by a straight line, 

{Ax'^{n)) = 2Dxjn. (13) 

The factor 7 is adopted from the definitions in [10]. Ex¬ 
pressions for the diffusion constant in the y-direction are 
obtained by replacing Xk by yk- Due to the slightly dif¬ 
ferent definitions, we call eqn. ( 12 ) short term diffusivity, 
and eqn. (13) the long term diffusivity. The averages are 
taken after initial transients have decayed. The number 
of periods one has to wait depends on the system size 
and the asymptotic behavior: While in the disordered 
regime, the system approaches its final state quickly, this 
time may be several thousand periods when the system 
approaches an ordered configuration. This shows up in 
5s, which can be used to distinguish ordered from disor¬ 
dered states. For disordered states, the diffusion constant 
is obtained by averaging over 10^ to 10^ periods. In the 
ordered states, the diffusion constant calculated from av¬ 
erages over shorter time intervals shows a drop towards 
zero. 

Since the simulation of large ensembles is computa¬ 
tionally expensive, we concentrate on a system size of 
N = 100 from here on. To verify that the results are not 
affected by the system size, we occasionally investigate 
systems of sizes N = 400 and 900 for a few parameter 


values. While the results for the diffusion constants seem 
to be independent of the system size as shown below, the 
time for transients to decay increases rapidly with system 
size and adds to the numerical challenges. 


C. Transitions for fixed Fi 

We begin the exploration of the parameter space of the 
system with a fixed interaction parameter T; = 0.1 and 
different shear rates Tg. We compare both short term 
and long term diffusivities in Fig. 6 . As anticipated from 
the results in the previous section, we notice two qualita¬ 
tively different states that can be distinguished by their 
diffusion constants. In one case, the system approaches 
a state with zero diffusivity and reversible motion. In 
the other case, the diffusivity does not vanish and the 
motion is irreversible. For this specific set of parame¬ 
ters, we computed the diffusivities both for a 100 and a 
900 particle system. Up to a shear rate of Tg just above 
30.0, the results are in very good agreement. Beyond 
this value, the system with 100 particles gives slightly 
higher values. This is most likely caused by finite size 
effects: At Tg = 36, which corresponds to a strain ampli¬ 
tude 7 o = 5.73, particles separated by Ay = •\/3/2 in the 
normal direction are advected by one box length relative 
to each other. Therefore, at the left and right turning 
points of the shearing motion, particles experience the 
same neighborhood. 

The transition between reversible and irreversible mo¬ 
tion as the shear rate is varied shows up rather clearly. 
For T; = 0.1, this transition takes place near Tg « 10.5. 
Fitting a power law of the form D = 5 (Fg — r 2 ,c)^ re¬ 
turns an exponent = 1.00, a critical value Ff^. = 9.6 
for the diffusivity in the shear direction, and /3j, = 0.47 
and ^ 2 c — 10-2 for the diffusivity in the normal direc¬ 
tion. The values for the critical shear rate are lower than 
the point where a non-zero diffusivity is first observed. 
The inset in Fig. 6 shows this critical region in more 
detail. A closer inspection of the diffusivities near the 
critical point reveals that they drop to zero abruptly. 
This suggests a small region where ordered and disor¬ 
dered motions coexist, so that the transition could be 
subcritical and hysteresis effects might be observed. Un¬ 
expected are the differences in the exponents: Parallel to 
the shear direction, the diffusivity grows linearly with the 
shear, but perpendicular to it, the diffusivity varies with 
a square root. Overall, the sharp transition is in quali¬ 
tative agreement with the results both from experiments 
and simulations, e.g. [ 10 , 11 , 16]. 

We find that the short term diffusivity eqn. (12) is 
slightly lower than its long term counterpart. Parallel to 
the shear direction, the difference approaches a constant 
value, whereas in the perpendicular direction, it dimin¬ 
ishes for larger shear rates. Our interpretation is that on 
short time scales, each particle is captured in a matrix of 
surrounding particles. Their presence introduces mem¬ 
ory and correlations into the process, thereby reducing 
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FIG. 6. (Color online) Diffusivities along the x- and y- 
direction, left and right panel respectively, for Fi = 0.1. Ad¬ 
ditional data of a 900-particle systems is marked by open 
symbols (short term: o; long term: o). For low shear rates, 
no diffusion is observed, corresponding to a reversible state. 
At a critical point 9.5 < r 2 ,c < 10.5 diffusions start to grow, 
both parallel and perpendicular to the shear, corresponding 
to an irreversible state. The inset in the left panel shows a 
magnification of the critical region. The short term diffusiv- 
ity is slightly smaller than the long term diffusivity. More¬ 
over, diffusion parallel to the shear is strongly enhanced and 
shows a different behavior: It grows linearly with increasing 
shear whereas diffusion in the perpendicular direction satu¬ 
rates. Boundary effects due to the system size seem to be 
negligible up to Fa « 30.0 since data points for 100 and 900 
particles coincide well. 

the diffusion coefficient. This effect is most pronounced 
in the parameter range close to the transition and be¬ 
comes weaker with increasing shear. 

D. Anisotropy in diffusion 

We also notice that horizontal diffusivities are much 
stronger than vertical ones. This is due to the shear act¬ 
ing only along the x-direction and was also observed in 
other investigations on similar systems [10, 16]. More¬ 
over, for larger shear rates, the vertical diffusivity seems 
to approach a fixed value whereas the horizontal dif¬ 
fusivity grows linearly over the investigated parameter 
range. This phenomenon of enhanced diffusion is known 
as advection-diffusion coupling and has been studied for 
the case of Brownian motion in a shear flow by Young 
et al. [34]. They showed that the diffusion constants are 
related to each other by 

Dx = Dx^o + (14) 

where 70 is the strain amplitude. Thus, parallel diffusion 
is enhanced by a coupling between normal diffusion and 
shearing. The bare longitudinal diffusivity without this 
effect is denoted Dx^. The mechanism can be illustrated 
in a three step process by decoupling diffusion and advec- 
tion: In the first step, we observe diffusion perpendicular 



FIG. 7. (Color online) Relation of the corrected parallel 
diffusivity (eqn. (14)) to the perpendicular diffusivity. The 
rescaled parallel diffusivity is still stronger, indicating a true 
anisotropy. Data was taken from the 900-particle system. 

to the shear so that particles are found in the layer above. 
In the next step, this layer is advected by the shear over 
some distance. In the last step, particles in the upper 
layer diffuse back to the original layer. Since they have 
been transported by the flow, they have traveled further 
than by horizontal diffusion alone. Equation (14) shows 
that for an isotropic system the bare longitudinal and 
transversal diffusivities coincide. 

In Fig. 7, we show the corrected parallel in relation to 
the perpendicular diffusivity. Data was taken with the 
900-particle system. We consider Dx/{1 + ^ 7 o), which 
in the limit of large 70 should approach Dy. We ob¬ 
serve that the bare parallel diffusivity is noticeably larger 
than the perpendicular diffusivity, contrary to the expec¬ 
tations for an isotropic system. Moreover, we can extract 
a factor between the two of approximately 1.5 over the 
whole range, even in the limit Tg « 50, i.e. | 7 g « 32. 
This suggests that the bare horizontal diffusivity Dxfi 
grows with shear rate. We confirmed that this behavior 
is not affected by boundary conditions by repeating the 
simulations in a quadratic box where we found the same 
increase. 

The differences between bare longitudinal and trans¬ 
verse diffusion disappear for stronger interactions Tj. The 
diffusivities for Tj = 1.0 in an ensemble of 900 parti¬ 
cles are presented in Fig. 8 . On the numerical side, we 
note that due to the strong interaction and the result¬ 
ing steeper gradients, integration times increase consid¬ 
erably. Below Tg = 28 the system has not yet converged 
towards an ordered state: While most particles are ar¬ 
ranged in a lattice configuration, unordered regions ex¬ 
ist which relax only very slowly. This is also reflected 
in the squared displacements: Instead of approaching 
a linear growth for larger 7 !, they obey a power law, 
(Aj/^) = 2Dy{'-ft)°‘ with a > 1. For larger shear rates, 
the simulation has converged and both the short term 
and long term diffusivity perpendicular to the strain co¬ 
incide. Virtually the same applies to the longitudinal 
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FIG. 8. (Color online) (Rescaled) diffusivities at Fi = 1.0. 
Below Fs = 28, the simulation has not converged towards the 
asymptotic ordered state: While most particles are arranged 
in a lattice structure, multiple defects exist which slowly grow 
out. Therefore, the long term diffusivity strongly deviates 
from the short term diffusivity and is neglected. At larger 
shear rates, we observe good agreement with eqn. (14), i.e. the 
anisotropy is caused by advection-diffusion coupling. 

diffusivity, though we do not observe perfect agreement. 
The cause for the deviation is the same as for Tg = 0.1, 
the particles being correlated for a few periods before 
showing diffusive behavior. The more important observa¬ 
tion is that the rescaled diffusivities D^.o in Fig. 8 nicely 
bracket the perpendicular diffusivity. Thus, for larger Tj 
the anisotropy can be completely explained by advection- 
diffusion coupling. 


E. Exploration of the Fi-Fs-parameter plane 

Next, we turn to the dependence on the shear rate 
and the interaction strength. In Fig. 9, both parallel and 
perpendicular diffusivities are shown in a large domain 
of the parameter space Fi x Fg. 

We observe a distinct boundary where diffusivities be¬ 
come finite as the shear rate Fg is increased. It nicely 
corresponds to the boundary determined by the loss of 
ordering in the system, observed in Fig. 4. For a bet¬ 
ter comparability, a smoothed representation of the curve 
where the ordering number jlkfcl drops below 0.4 is shown 
in Fig. 9 as well. 

In general, the diffusivity increases with the interaction 
parameter Fj up to a critical value where it drops to zero. 
A small diffusivity D for small values of the interaction 
parameter F; is plausible, since a relatively higher shear 
rate Fg is needed to bring particles closer together and 
to allow the system to display large spatial fluctuations. 
The dominant force is the external forcing, and after one 
period the particles are very close to where they were one 
period before. This leads to a small diffusion constant. 
In fact, the relaxation time becomes so long that we can¬ 
not measure a proper diffusion constant. However, there 


FIG. 9. (Golor online) The long term diffusion constants 
Dx (left panel) and Dy (right panel) for different values of 
Fi and Fg. The diffusivity is obtained by fitting data with 
eqn. (13). The black solid lines indicate Dx = 1 x 10“® 
and Dy = 1 x 10“®, respectively. The dashed lines are the 
smoothed curves taken from Fig. 4, indicating the loss of 
structural ordering. Results for Fg > 5 • 27r have to be treated 
with care due to the limited box-width of Lx = 10: At this 
point, during one period each particle slips beyond five par¬ 
ticles in each direction. This leads to a stabilizing effect for 
Fi > 0.5: Although the shear force is strong, at the turning 
points of the oscillation the particle interaction still dominates 
and promotes ordering into a lattice configuration. This be¬ 
havior is not observed for larger box-sizes. 


are qualitative differences depending on the shear rate 
Fg: For small shear rates, the particles are evenly dis¬ 
tributed, but unordered, whereas for larger shear rates 
the particles are uniformly distributed. This is also sup¬ 
ported by two-particle correlation functions. However, 
the transition between the two regimes is not as sharp as 
for stronger particle interactions. For large interaction 
parameters we see a similar behavior, but for a different 
reason. Since the interaction is comparatively strong, we 
need a high shear rate to keep the particles away from 
their equilibrium positions. On the other hand, for large 
Fj, the settling time scale is shorter, and since the period 
of the external forcing is fixed to 1 , we expect to find a 
critical region where particles settle into an ordered state 
before they can start moving in a random fashion. This 
is reflected in a relatively sharper transition in Fig. 9 for 
Fi « 1 . As would be anticipated, stronger shear rates are 
required as the particle interaction is increased. 

For a dilute solution of particles, Pine et al. [10] found 
that the diffusivity is independent of the period of the 
forcing. In order to check whether our results are in¬ 
dependent of the driving period, we have to vary the 
period while keeping the strain 7 fixed. With our dimen¬ 
sionless parameters, this corresponds to a variation of F;, 
while Fg is kept fixed. From Fig. 9, we see that for small 
Fi < 0 . 2 , i.e. short enough periods and weak interactions, 
the onset of a finite diffusivity and the loss of ordering 
is independent of the driving period. For longer periods 
and stronger forces though, the interaction between par- 























tides has more time to drive the system to the ordered, 
reversible state. 

However, the relation to the experiments of Pine et al. 
[10] is more complicated. In the experiments, the shear 
rate was fixed and the shear amplitude was changed by 
varying the period. This implies a change in Pi as well, 
because the particles have more time to interact. In the 
fluids context, the interaction comes from the hydrody¬ 
namic forces, and they become smaller if the shear rate 
(and thus the speed) is reduced. The parameter regime 
where both systems show similar phenomenology then is 
the range of small Pi, where the transition is determined 
by Pg alone. 



IV. LINEAR STABILITY ANALYSIS 

Investigation of the long term trajectories of simula¬ 
tions in the reversible regime show that they approach 
regular lattices. That regular lattices are force equilib¬ 
ria, even under periodic shear, follows from the fact that 
to each particle pair there is a corresponding one with 
opposite forces. More formally, consider a lattice config¬ 
uration with a point symmetry in particle spacings, i.e. 
a situation where for every pair i, j one can find another 
pair i, k, such that Then it is easy to see 

that the inter-particle forces in eqn. (6) vanish. Since the 
shear force only introduces an affine transformation, the 
force balance stays unaffected during the shear process. 
Therefore, the regular lattices under shear correspond to 
a periodic cycle of system (6), 

x°(t) = sin( 27 rt)e„ (15) 

for all parameter values Pi and Pg. The initial regular lat¬ 
tice points are denoted and their time-dependent 

cousins are x°(t). The regular lattices we consider are 
the ones obtained under small shear (Fig. 10): the equi¬ 
lateral triangular lattice with particles aligned along the 
shear direction, the same lattice rotated by 90°, and the 
rectangular lattice aligned with the shear. In order to 
determine their stability against small perturbations, we 
have to perform a Floquet analysis, on account of the 
periodic variation of the sheared lattices. 


A. Floquet analysis 

We linearize the evolution equation (6) around the pe¬ 
riodic cycle (15) 


x,{t) = x°{t) + 6i{t). (16) 

This yields a system of linear ordinary differential equa¬ 
tions for the perturbation vector d{t), 

^ = Fi {V{t) - M{t)) 8 + r,S{t)5 = W{t)S. (17) 
at 


FIG. 10. (Color online) The three particle configurations 
considered in the linear stability analysis, from left to right: 
equilateral triangular lattice aligned along the x- and j/-axis 
and rectangular lattice. These are also the asymptotic states 
found in the dynamical system. 


The right-hand side of this system consists of a time- 
periodic coefficient matrix yV’(t) composed of 2N x 2N 
matrices 


and 


m = 




M{t) = 



-Si,2 ■ 


Bi^2 

0 • 

• B2,N 


B2,N ■ 

■ 0 / 


(18) 


(19) 


where Gi{t) = ^i,j (^) with the 2x2 Jacobi matrices 
of the pairwise particle interactions. 





4(x° 



( 20 ) 


Here, 0 denotes the dyadic product and I 2 is the 
identity-matrix in two dimensions. Finally, we have the 
Jacobi matrix of the external forcing. 



■) 

SM = (g 

V 0 


\ / 


( 21 ) 


Since the coefficient matrix has a periodic time- 
dependence, the linear stability analysis requires an in¬ 
tegration of the equations over a full period [35]. To 
this end, we compute the principal fundamental matrix 
$(t) of system (17) with initial condition <i)(0) = ll 2 Ar. 
The matrix after a full period, 4)(1) = $, is the Floquet 
matrix and is not necessarily a symmetric matrix. The 
stability properties of system (17) then depend on the 
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eigenvalues of $ only. The eigenvalues cr of $ are called 
Floquet multipliers. For stability, we require |cr| for all 
possible eigenvalues cr of $ to be less than one. There are 
two neutral modes (tr = 1), connected with the invari¬ 
ance of perturbations along the orbit and translational 
invariance in the x-direction. 

As an important technical detail we mention that the 
boundary conditions have to be implemented such that 
they are compatible with the symmetries of the system 
so that the symmetries of the matrices are preserved. To 
this end, we introduce an interaction radius around each 
particle, and take all interactions with particles inside the 
circle into account. The radius is chosen to be at least of 
the order of the box dimensions. 


B. Stable and unstable configurations 

We studied the stability of configurations with N = 
100,400 and 900 particles. We find that results are qual¬ 
itatively the same, with the stable modes being closer to 
unity for more particles. We investigated three different 
configurations as described above, and for the rotated 
hexagonal configuration we swapped and Ly. 

It turns out that for F, = 0.1, at least one configuration 
is stable for values of Fg < 10.5. This value corresponds 
to a shear rate where each particle is shifted by 1.5 par¬ 
ticle distances in relation to the neighboring rows. For 
larger values, the periodic cycles remain unstable. The 
largest multipliers are shown in Fig. 11 for a wide range of 
the shear rate. We find three regions with different stable 
configurations. For Fg < 3.5, both the hexagonal state 
and its rotated counterpart are stable. Thereafter, the 
rectangular state is stable up to a shear rate of Fg < 7.5 
until finally the hexagonal configuration is stable again. 
The changes occur at shear rates corresponding to rela¬ 
tive displacements of neighboring rows of approximately 
0.5, 1 and 1.5 particle distances. They nicely correspond 
to the modulations of the orientation number shown in 
Fig. 4. Only in a small interval at Fg « 7.5, all consid¬ 
ered lattices appear to be unstable. The inspection of 
long term trajectories shows, that the system relaxes to¬ 
wards a mixture of both the hexagonal and rectangular 
lattice and still is in the reversible regime. 

The different behavior of the rotated hexagonal con¬ 
figuration can be understood from a geometrical point 
of view: In the hexagonal lattice with particles aligned 
along the shear direction, the initial configuration is re¬ 
produced after a shear displacement by one particle dis¬ 
tance. In the rotated lattice a shift by four particle dis¬ 
tances is needed before the lattice returns, corresponding 
to a shear rate of Fg « 20. An analysis of the Jacobian 
indicates that it picks up strongly unstable modes during 
this process. 

We do not find any qualitative differences between the 
400 and 900 particle ensemble. Quantitatively, the stable 
modes become less stable with increasing number of par¬ 
ticles. Additionally, the points where the configurations 


(a) 



Fs 




FIG. 11. (Color online) The largest eigenvalues Umax of $ 
in dependence on Fg for the three investigated configurations 
with N = 100, 400, and 900 particles from top to bottom. 
Neutral modes are omitted. In each upper panel, markers 
indicate eigenvalues |crniax| > 1 and thus unstable configu¬ 
rations, and stable modes (Irrmaxl < 1) in the lower panels 
accordingly. In between, the regions of stability are repre¬ 
sented by colored bars. We find that for values of Fg < 11, 
linearly stable states exist. Beyond, linear stability is lost. 
This is in good agreement with the diffusion constant calcu¬ 
lated earlier. With increased number of particles, the stable 
modes become less stable, the qualitative difference between 
400 and 900 particles being rather marginal. The stable patch 
at Fg > 12 in the 100-particle system can clearly be identified 
as a finite size effect due to the limited box-width. The value 
of Fi is set to 0.1 in all calculations. 
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change stability are slightly shifted. The magnitudes of 
the unstable modes, however, are almost identical. The 
100 particle system, on the other hand, clearly shows 
signs of boundary effects. In the range 12 < Tg < 14, the 
rectangular lattice becomes stable again. This is neither 
observed in the larger systems, nor in the simulations. 
The source might be the limited box width in conjunction 
with the long range potential. Except for this detail, the 
100 particle system qualitatively agrees with the larger 
ensembles as well. 

A variation of the interaction strength r, reveals that 
the location of stable regions is independent thereof for 
Ti < 0.2. For a larger interaction strength, more stable 
regions emerge, in general accordance with Fig. 9. How¬ 
ever, the computational demand considerably increases 
as Fi is increased, so that a systematic investigation with 
an adequately large number of particles is unfeasible. 


V. PARTICLE DISTRIBUTIONS 

A. Radial distribution functions 


Additional information about the structure of the sys¬ 
tem that goes beyond the value of the diffusion constants 
and reversibility can be obtained through correlation 
functions that contain information about the arrange¬ 
ments of particles. The radially averaged two-particle 
correlation function counts the number of particles within 
a ring bounded by the radii r — 6r and r Sr, normalized 
by the mean number of particles in such a ring. 


_ 1 Ir-Sr/2 S{r,j - f)dr 

^ Ir-Sr/2 


( 22 ) 


where p = N/{LxLy) is the average particle density of the 
system. For a crystalline structure, one expects sharp 
peaks corresponding to the underlying lattice, the first 
one at the next neighbor distance, and so on. In a liquid 
state no long range ordering is found so that the correla¬ 
tion approaches one at larger distances. As neighboring 
particles repel each other, one expects a gap close to the 
particle and a strong peak at the position of the nearest 
neighbors. Fig. 12 shows the correlation functions of the 
system for different shear rates Fg (and F = 0.1). To 
obtain better statistics all simulations in this section are 
done for systems of 900 particles. Inspection of the final 
states and the displacements per particle show that they 
behave exactly as the smaller systems of 100 particles. 
Each ensemble is taken at a full period after a total sim¬ 
ulation time of 11000 periods, i.e. we assume the system 
has reached its asymptotic state, which is certainly the 
case for the systems at larger shear rates. 

As expected, at low shear rates the ordered structure 
of the state is reflected in the correlation function. At 
Fg = 6.0, we hnd two superimposed peaks at r = •\/3/2 
and 1.0, corresponding to the rectangular lattice. At a 



r 


FIG. 12. (Color online) Radial correlation functions of 
ensembles at different shear rates. For clearer distinction, 
graphs are shifted by a constant. At low shear rates (Fg = 6.0 
and 8.0, respectively), the crystalline structure is reflected by 
isolated peaks at the next neighbor distances; At Fg = 6.0 the 
rectangular lattice with a — 1.0 and b = \/3/2, and the hexag¬ 
onal lattice with a = 1.0 at Fg = 8.0. At higher shear rates, 
in the irreversible regime, the correlation function shows fea¬ 
tures of a liquid: Beyond a minimum distance, the function 
quickly approaches 1.0 and displays only small fluctuations. 
Simulation of 900 particles, Fi = 0.1. 


slightly increased shear rate of Fg = 8.0 the asymptotic 
state changes, as can be seen in the correlation function 
which exhibits isolated peaks at r = 1.0 and y/i, corre¬ 
sponding to the next neighbor distances in a hexagonal 
lattice. At higher shear rates after the motion became ir¬ 
reversible, the correlation function shows features of a liq¬ 
uid: Beyond a minimum distance owing to the repelling 
potential, the function quickly approaches 1.0 and shows 
only small modulations. For even larger shear rates, the 
correlation function becomes featureless except for the 
next neighbor repulsion: the system is in a gas-like state. 
However, no spatial information can be extracted in the 
latter cases. The next neighbor distance is identified as 
0.8 at Fg = 16.0, and it decreases to 0.7 at Fg = 30.0. 
This points to an improvement of the mixing process with 
increasing shear rate, which leaves the particles less time 
to relax and to return to larger separations. 


B. 2-d correlation functions 

In order to obtain more spatial information, we now 
turn to the spatially resolved two-particle correlation 
function 


y+^y/"^ 

g{x,y)=^ (23) 

” v-Sv/2 b 


where n = NpSxSy. Just as its radial counterpart, the 
2d correlation function relates the number of particles in 
a rectangular box of size 6x x 6y around the distance x 
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FIG. 13. (Color online) The 2-dimensional correlation func¬ 
tion g{x, y) at low shear rates. Colors indicate the neighbor 
density. At the white patches, g{x, y) exceeds 2.0. Just as 
the radial correlation function g{r), g{x,y) reflects the crys¬ 
talline structure: At Fg = 6.0, (a), we observe the rectangular 
lattice structure, whereas at Fg = 8.0, (b), the hexagonal lat¬ 
tice is recovered. In both cases, the blur is caused by minor 
lattice defects whereby in the latter case almost perfect ver¬ 
tical alignment is observed. The white triangles indicate how 
far the populated rows have been sheared against each other 
over the past quarter period. They are set at arbitrary height 
such that features of the correlation function are not occluded. 
Simulation of 900 particles, Fi = 0.1 


X 



-4 -2 0 2 4 

X 


and y to the mean number of particles expected in such 
a region. 

Figure 13 shows the 2d correlations for the ordered 
configurations at low shear rates. Data was acquired in 
the same manner as for the radial correlation function. 
In the plot, the correlation function is color coded from 
dark brown to white with increasing density and cen¬ 
tered around the neutral density of unity. As expected 
from the radial correlation, for Fg = 6.0 we find a rectan¬ 
gular lattice with most particles scattered closely to the 
lattice points and a few particles found on the grid lines. 
This is due to some minor lattice defects. Accordingly, 
at F = 8.0 we recover the hexagonal lattice. Again, par¬ 
ticles gather closely to the lattice sites, at which vertical 
alignment is almost perfect. The horizontal scatter is due 
to over- or under-population in rows which locally lead 
to rectangular arrangements of particles. Note that both 
ensembles show mirror symmetries along the x- and y- 
axis, i.e., there is no memory in the system if it has been 
sheared to the left or to the right over the last half period. 

The picture changes upon advancing to higher shear 
rates. In Fig. 14 are shown the 2d correlations g{x,y) 
across half a period at Fg = 16.0 and 30.0, respectively. 
In figure (a) are shown snapshots of Fg = 16.0 at the 
beginning of a period, at the turning point of f = T/4 
when the flow has stopped and strain is maximal, and af¬ 
ter half a period at t = T/2, when the flow has reversed 
its direction and strain is zero again. The most promi¬ 
nent feature is the break-down of the symmetry. At the 
start of the period, t = 0, we observe a dark brown rhom¬ 
bic region around the origin where no other particles are 
found. It is oriented to the right. As would be expected 


FIG. 14. (Color online) The 2-dimensional correlation func¬ 
tion g{x, y) at high shear rates over half a period at (a) 
Fg = 16.0 and (b) Fg = 30.0. Colors indicate the particle 
density. At the white patches, g(x,y) exceeds 2.0. Each sub¬ 
figure shows, from top to bottom, snapshots at t = 0, t = T/4, 
and T = T/2. We observe a break-down of the mirror symme¬ 
try found at lower shear rates. While the free space is tilted to 
the right at the beginning of a period, half a period later it is 
tilted to the opposite direction. Moreover, we Hnd remnants 
of the hexagonal lattice in the form of stripes of increased den¬ 
sity spaced by roughly the layer distance of the crystal. The 
white triangles indicate how far the populated rows have been 
sheared against each other over the past quarter period. They 
are set at arbitrary height such that features of the correla¬ 
tion function are not occluded. Simulation of 900 particles, 
Fi = 0.1 


from the symmetry of the flow, half a period later we 
observe a similar rhombus, but this time oriented to the 
left. In between, at the turning point of the flow, its 
outer boundary at y ~ ±0.5 is stretched along the strain 
even more, while at the centerline j/ = 0 it keeps its 
boundaries. In effect, this leads to less populated protru¬ 
sions above and below the centerline, concurrent with the 
relative strain. They closely resemble the gaps found in 
models describing the emergence of moonlet propellers in 
Saturn’s A ring [36]. A ‘propeller-shaped’ region was also 
described by Keim et al. [37] in their study of the model 
in [18]. Due to a different definition of the shear protocol, 
our snapshot at t = T/4 corresponds to f = T/2 in their 
system, and by symmetry relates to t = T. A remarkable 
difference to their observation is that they describe the 
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FIG. 15. Formation of a shear rate dependent depletion zone (a) and ratio of the bare diffusivities and length scales in 
dependence on the shear rate (b). In (a), we consider two particles (•), one being fixed at the origin, and the other one being 
subject to the forces (arrows) due to particle interaction Fij at Fi = 0.1 and the shear Fghear, which linearly increases with the 
separation in the t/-direction. Black lines represent the points where horizontal forces are balanced. Below, the free particle 
moves to the left, i.e. particles outside the area cannot enter it. With increasing shear rate, the maximum value i/max decreases. 
The red line marks the position of largest thickness of the depletion zone, changing with Fg. In (b), the ratio of the bare 
diffusivities Dy/Dj;^ shows the same behavior as the quotient j/max/Fg. 


correlations for the reversible state, whereas we here find 
it in the irreversible state. Therefore, also the irregular 
state maintains a phase memory of the shear process, at 
least for the interactions studied here. 

We can understand this phenomenon when starting at 
t = —TjA. At this time, the flow is at its left turning 
point. We now look at the line a; = 0. At short distances, 
particles on the left of this line are repelled to the left, 
and over the next quarter period, the flow is not fast 
(i.e. strong) enough to pull them beyond a: = 0. Thus, 
at t = 0 we observe the rhombus: it simply corresponds 
to the line a: = 0, advected over the last quarter period. 
This suppression of transport continues up to t = T/4. 
Since this does not apply to the centerline, the rhombic 
shape gets distorted and we observe the aforementioned 
protrusions. As soon as the low density area leaves the 
core area, particles from the other sides are pushed into 
the void. This is reflected by a slight increase in g{x,y) 
towards the tips of the protrusions. However, at t = T/2, 
the remnants are still visible (brown (gray) areas in the 
upper right and lower left). 

The maxima of the correlation function are found at 
X « ±1.0 and ?/ = 0 at every time step. This corresponds 
to a strong localization of next neighbors parallel to the 
shear. Additionally, in the perpendicular direction we 
find a modulation in particle density, with maxima at 
y Ki 0 and ±0.6 ... 0.8, and multiples thereof. The mod¬ 
ulations extent over several particle distances parallel to 
the shear before they become weaker and decay. This has 
two implications. First, although the system is far from 
its equilibrium of a hexagonal lattice configuration, it still 
shows the vertical density modulation corresponding to 
the layers in the lattice. With increasing shear rate, this 
becomes less pronounced, but is still detectable. Second, 
those maxima in the correlation function hint at a gen¬ 
eral anisotropy: The next neighbor distance parallel to 


the shear is conspicuously larger than the next neighbor 
distance in the perpendicular direction, as reflected in 
the location of the layers. 

At the higher shear rate Tg = 30.0, shown in Fig. 14(b), 
the picture qualitatively stays the same. We still observe 
the density modulations and the flow dependent varia¬ 
tion of the next neighbor distances. What changes is the 
shape of the particle-free space. At the beginning of the 
period, it gets stretched even more, but does not resemble 
a rhombus any more. This is due to the fixed neighbor 
relations along the centerline at x « ±0.8... 1.2. At the 
turning point, t = T/4, the free space is a slightly tilted 
oval instead of a rounded rhombus. As expected, the 
density modulations in the perpendicular direction have 
a smaller wavelength than at Fg = 16.0, and the first 
maxima are found at j/ « ±0.5 ... 0.6. 

The reduction of the vertical spacing with shear rate 
can be explained within a simple model. To this end, we 
consider two particles, a fixed one placed at the origin 
and a free one, as sketched in Fig. 15(a). For now, we re¬ 
place the periodic shear by a constant one which linearly 
increases with the separation y. Thus, the free particle is 
subject to two different forces, the inter-particle force 
and the shear force Fghear acting along the x-direction. 
Taking only forces parallel to the shear into account, we 
can compute the particle position where horizontal forces 
balance, 

Tghear = — Tij,||. (24) 

If the mobile particle is close to the x-axis, forces are 
balanced at large distances. For larger separations in y, 
the shear force dominates and the force from the station¬ 
ary particle is too weak to balance them. The curves 
of horizontal force equilibrium are shown in Fig. 15(a) 
for various shear rates. Particles in the area below the 
curves move to the left and are deflected upwards un- 
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til they are finally far enough away to pass the station¬ 
ary particle. This deflection can also be observed in the 
correlation function in Fig. 14 at t = T/d, where the 
maximum at a; = ±1 is bent away from the centerline. 
Due to the periodic repetition and the many-particle in¬ 
teractions, this eventually leads to the formation of a 
somewhat larger depletion zone where g{x, y) is very low. 
When increasing the shear rate, forces are balanced closer 
to the origin and so the area of the depletion zone de¬ 
creases. The maximum separation ymax varies along a 
straight line with Omax = 27r/3 the angle towards the x- 
axis, and approaches the origin in the limit of very large 
shear rates. By contrast, the hard-core interactions used 
by Keim et al. [37] do not lead to a shear rate depen¬ 
dent separation: They report a constant separation in 
the perpendicular direction with boundaries at y = ±1. 

Advection-diffusion coupling (14) is based upon the as¬ 
sumption that a single particle performs diffusive motion. 
However, in the system studied in this paper we have a 
multi-particle interaction and thus the effects of the sur¬ 
rounding particles have to be taken into account. These 
interactions determine the anisotropies of the diffusivities 
in Fig. 7. To see this, note that the free path of a par¬ 
ticle is determined by the regions of low particle density, 
i.e. low g{x,y). From the previous observations and the 
simple model, we can derive two assumptions: First, the 
free path perpendicular to the shear is governed by the 
maximum separation j/max- Obviously, they do not cor¬ 
respond directly to each other since j/max is smaller than 
the separation in Fig. 14. Yet, ymax should determine 
the general dependency on the shear rate. Second, the 
free path parallel to the shear is set by the low-density 
protrusions, which roughly increase along with the shear 
rate Fg. 

These assumptions then suggest that the ratio of the 
bare diffusivities Dy/D^fi should be proportional to 
ymax/bs. Fig. 15(b) shows that this is indeed the case 
for Fg > 15. This shows that the unexpected anisotropy 
of the system is connected with the emergence of two 


different length scales which change with the shear rate. 


VI. SUMMARY 

We have seen that even a model as simple as (6) dis¬ 
plays features usually associated with much more com¬ 
plex systems. We found chaotic and ordered behavior 
which can occur in different regions of the parameter 
space (ri,rg): The transition is similar for each value 
of the interaction strength: For small shear strains, after 
some transient time the system becomes reversible, the 
shear can be considered as a small perturbation. Beyond 
an interaction-dependent critical shear rate, this feature 
is lost and we observe chaotic motion. For stronger in¬ 
teractions, the critical point is delayed to higher shear 
rates. 

We identified the diffusivity as an indicator of chaotic 
motion. As would be expected for a diffusive system 
subject to shear, we observe advection-diffusion cou¬ 
pling: The parallel diffusivity is strongly enhanced as 
compared to perpendicular diffusivity. For large inter¬ 
action strengths, we can explain this anisotropy using 
eqn. (14). For low interactions though we observe addi¬ 
tional anisotropic effects, which we can attribute to the 
2d-correlations of particles: Whereas the free space par¬ 
allel to the shear increases with the shear rate, the next 
neighbor distance perpendicular to the shear decreases 
with increasing shear rate. Both effects combined allow 
us to explain the relative anisotropy of the bare diffu- 
sivities. Additionally, we were able to identify a phase 
dependency of the particle distributions, which is con¬ 
cealed when looking upon e.g. stroboscopic maps. 

Finally, we could relate the critical shear rate associ¬ 
ated with the loss of reversibility to the stability proper¬ 
ties of sheared regular lattices. So at least for the present 
system the onset of irreversibility is connected with an 
instability of the asymptotic, regular state. 


[1] E. L. Hahn, Phys. Rev. 80, 580 (1950). 

[2] H. Y. Carr and E. M. Purcell, Phys. Rev. 94, 630 (1954). 

[3] G. M. Homsy, ed., Multimedia Fluid Dynamics, 2nd ed. 
(Cambridge University Press, 2008). 

[4] W. Niggemeier, G. von Plessen, S. Sauter, and 
P. Thomas, Phys. Rev. Lett. 71, 770 (1993). 

[5] J. Chaiken, R. Chevray, M. Tabor, and Q. M. Tan, Proc. 
R. Soc. A 408, 165 (1986). 

[6] J. A. G. Roberts and G. R. W. Quispel, Phys. Rep. 216, 
63 (1992). 

[7] G. Casati, B. V. Chirikov, I. Guarneri, and D. L. Shep- 
elyansky, Phys. Rev. Lett. 56, 2437 (1986). 

[8] H. M. Pastawski, P. R. Levstein, and G. Usaj, Phys. 
Rev. Lett. 75, 4310 (1995). 

[9] B. Eckhardt, J. Phys. A 36, 371 (2003). 

[10] D. J. Pine, J. P. Gollub, J. F. Brady, and A. M. Leshan- 
sky. Nature 438, 997 (2005). 


[11] J. S. Guasto, A. S. Ross, and J. P. Gollub, Phys. Rev. 
E 81, 061401 (2010). 

[12] B. Metzger and J. E. Butler, Phys. Fluids 24, 021703 

( 2012 ). 

[13] R. Jeanneret and D. Bartolo, Nat. Commun. 5, 3474 
(2014). 

[14] A. Franceschini, E. Filippidi, E. Guazzelli, and D. J. 
Pine, Phys. Rev. Lett. 107, 250603 (2011); Soft Matter 
10, 6722 (2014). 

[15] N. C. Keim and P. E. Arratia, Soft Matter 9, 6222 (2013); 
Phys. Rev. Lett. 112, 028302 (2014). 

[16] B. Metzger and J. E. Butler, Phys. Rev. E 82, 051406 

( 2010 ). 

[17] G. During, D. Bartolo, and J. Kurchan, Phys. Rev. E 
79, 030101 (2009). 

[18] L. Corte, P. M. Chaikin, J. P. Gollub, and D. J. Pine, 
Nat. Phys. 4, 420 (2008). 



14 


[19] L. Corte, S. J. Gerbode, W. Man, and D. J. Pine, Phys. 
Rev. Lett. 103, 248301 (2009). 

[20] G. I. Menon and S. Ramaswamy, Phys. Rev. E 79, 061108 
(2009). 

[21] N. C. Keim and S. R. Nagel, Phys. Rev. Lett. 107, 010603 

( 2011 ). 

[22] L. Mohan, C. Pellet, M. Cloitre, and R. Bonnecaze, J. 
Rheol. 57, 1023 (2013). 

[23] 1. Regev, T. Lookman, and C. Reichhardt, Phys. Rev. 
E 88, 062401 (2013). 

[24] D. Fiocco, G. Fofii, and S. Sastry, Phys. Rev. E 88, 
020301 (2013); Phys. Rev. Lett. 112, 025702 (2014). 

[25] S. Slotterback, M. Mailman, K. Ronaszegi, M. van Hecke, 
M. Girvan, and W. Losert, Phys. Rev. E 85, 021309 
( 2012 ). 

[26] G. F. Schreck, R. S. Hoy, M. D. Shattuck, and C. S. 
OHern, Phys. Rev. E 88, 052205 (2013). 

[27] J. R. Royer and P. M. Ghaikin, Proc. Natl. Acad. Sci. 
112, 49 (2014). 

[28] N. Mangan, C. Reichhardt, and G. J. Olson Reichhardt, 
Phys. Rev. Lett. 100, 187002 (2008). 

[29] W. Zhang, W. Zhou, and M. Luo, Phys. Lett. A 374, 
3666 (2010). 

[30] S. Okuma, Y. Suzuki, and Y. Tsugawa, Phys. C 470, 


S842 (2010); A. Motohashi and S. Okuma, J. Phys. Gonf. 
Ser. 302, 012029 (2011); S. Okuma, Y. Tsugawa, and 
A. Motohashi, Phys. Rev. B 83, 012503 (2011). 

[31] A. Melzer, A. Homann, and A. Piel, Phys. Rev. E 
53, 2757 (1996); A. Melzer, V. A. Schweigert, 1. V. 
Schweigert, A. Homann, S. Peters, and A. Piel, Phys. 
Rev. E 54, R46 (1996); V. A. Schweigert, 1. V. 
Schweigert, A. Melzer, A. Homann, and A. Piel, Phys. 
Rev. Lett. 80, 5345 (1998). 

[32] A. W. Lees and S. F. Edwards, J. Phys. C: Solid State 
Phys. 5, 1921 (1972). 

[33] See Supplemental Material at URL for movies. 

[34] W. R. Young, P. B. Rhines, and C. J. R. Garrett, J. 
Phys. Oceanogr. 12, 515 (1982). 

[35] C. Chicone, Ordinary Differential Equations with Appli¬ 
cations, Texts in Applied Mathematics, Vol. 34 (Springer 
New York, 2006). 

[36] F. Spahn and M. Sremcevic, Astron. Astrophys. 358, 368 
(2000); M. Sremcevic, F. Spahn, and W. J. Duschl, Mon. 
Not. R. Astron. Soc. 337, 1139 (2002); M. Sremcevic, 
J. Schmidt, H. Salo, M. Seiss, F. Spahn, and N. Albers, 
Nature 449, 1019 (2007). 

[37] N. C. Keim, J. D. Paulsen, and S. R. Nagel, Phys. Rev. 
E 88, 032306 (2013). 


